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Abstract 

Background: The ability to perform quantitative studies using isotope tracers and metabolic flux analysis (MFA) is 
critical for detecting pathway bottlenecks and elucidating network regulation in biological systems, especially those 
that have been engineered to alter their native metabolic capacities. Mathematically, MFA models are traditionally 
formulated using separate state variables for reaction fluxes and isotopomer abundances. Analysis of isotope 
labeling experiments using this set of variables results in a non-convex optimization problem that suffers from 
both implementation complexity and convergence problems. 

Results: This article addresses the mathematical and computational formulation of 13 C MFA models using a new 
set of variables referred to as fluxomers. These composite variables combine both fluxes and isotopomer 
abundances, which results in a simply-posed formulation and an improved error model that is insensitive to 
isotopomer measurement normalization. A powerful fluxomer iterative algorithm (FIA) is developed and applied to 
solve the MFA optimization problem. For moderate-sized networks, the algorithm is shown to outperform the 
commonly used 13CFLUX cumomer-based algorithm and the more recently introduced OpenFLUX software that 
relies upon an elementary metabolite unit (EMU) network decomposition, both in terms of convergence time and 
output variability. 

Conclusions: Substantial improvements in convergence time and statistical quality of results can be achieved by 
applying fluxomer variables and the FIA algorithm to compute best-fit solutions to MFA models. We expect that 
the fluxomer formulation will provide a more suitable basis for future algorithms that analyze very large scale 
networks and design optimal isotope labeling experiments. 



Background 

Metabolic Pathway Analysis 

Metabolism is the complete set of chemical reactions tak- 
ing place in living cells. These chemical processes form 
the basis of all life, allowing cells to grow, reproduce, 
maintain their structure and respond to environmental 
changes. Metabolic reactions are divided into groups 
called metabolic pathways, which are typically con- 
structed heuristically according to their connectivity and 
presumed function [1]. Each metabolic pathway is char- 
acterized by a set of chemical reactions that transform 
substrates into end products while generating intermedi- 
ate byproducts. Due to its importance in medicine and 
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biotechnology, metabolic pathway research has become a 
highly active field of investigation [2]. 

Initially, the structure of metabolic pathways was 
examined by identifying their intermediate compounds. 
Subsequently, the various biochemical reactions con- 
necting these compounds were mapped. Due to the suc- 
cess of this research, the topological structure of many 
metabolic pathways is nowadays fully documented [3]. 
The next step in the progression of metabolic pathway 
research involves quantification of the rates of these var- 
ious chemical reactions, also known as "fluxes". The 
values of these rates are affected by various environmen- 
tal conditions and can change rapidly in response to 
perturbations. Nevertheless, if the environmental para- 
meters are held fixed and stable, the network can attain 
a steady state in which the concentrations of all network 
metabolites are assumed constant over time. This, of 
course, implies that the rates of their input and output 
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reactions must balance. The latter imposes a set of lin- 
ear constraints on the metabolic fluxes, known as "stoi- 
chiometric balance equations" [4]. Unfortunately, since 
the number of unknown fluxes typically exceeds the 
number of independent stoichiometric balances, these 
constraints are insufficient to completely identify the 
metabolic network. In order to overcome this lack of 
information, additional constraints must be provided to 
the stoichiometric mathematical model to estimate the 
values of the network fluxes [5]. 

13 C Isotope Labeling Experiments 

Various experimental techniques have been developed to 
enable measurement of intracellular metabolic fluxes, 
either directly or indirectly. One of these approaches 
makes use of isotope labeling experiments. In this 
method, the metabolic system is administered a known 
amount of an isotopically labeled substrate (such as glu- 
cose labeled with 13 C at specific atom positions). By 
measuring the resulting labeling patterns of intracellular 
metabolites after steady state has been achieved, addi- 
tional flux information is obtained. 

One major drawback of this experimental approach is 
the high complexity and computational intensity of the 
metabolic flux analysis (MFA) required to interpret 
these labeling measurements. In their series of articles, 
Wiechert et al. [6-9] constructed a systematic approach 
for performing this analysis. They show that measure- 
ments of the relative abundance of various isotope iso- 
mers, also known as "isotopomers", contain enough 
information to fully identify the metabolic fluxes of the 
network. Formulating the problem using isotopomer 
variables (or a transformed set of isotopomer variables 
referred to as "cumomers"), Wiechert et al. posed the 
flux estimation problem as a non-convex least-squares 
minimization, assuming random error is added to their 
isotopomer measurements. The resulting high-dimen- 
sional non-convex problem suffers from various draw- 
backs, such as slow convergence and notable probability 
of attaining local minima. Several optimization algo- 
rithms have been developed in order to address these 
drawbacks. Early approaches used iterative parameter- 
fitting algorithms [8], evolutionary algorithms [10] and 
simulated annealing [11]. Furthermore, several investiga- 
tions have been conducted in order to assess the accu- 
racy of these results [9,12,13]. Recently, a novel method 
to decompose the metabolic network into Elementary 
Metabolite Units (EMUs) was introduced [14] and 
implemented into the OpenFLUX software package [15]. 
This decomposition effectively reduces the size of the 
optimization problem by efficiently simulating only 
those isotopomers that contribute to the measurement 
residuals. Nevertheless, all of these algorithms suffer 
from several major drawbacks due to the standard 



isotopomer-flux variables used in formulating the opti- 
mization problem: 

♦ Presence of unstable local minima: due to the non- 
convex nature of the objective function. 

♦ Complex mathematical representation and compu- 
tational implementation. This results in the need for 
ad-hoc algorithms and mathematical analysis, and 
long running times are required for reliable 
convergence. 

The OpenFLUX implementation, for example, may 
require several dozens of convergence iterations with 
various initial values in order to achieve acceptable 
probability of obtaining the optimal set of fluxes in any 
one of its attempts. In addition, due to the chosen 
objective function, it is also commonly required to esti- 
mate scaling factors for each isotopomer measurement, 
because of the fact that the available experimental tech- 
niques are only capable of measuring isotopomer frac- 
tions up to a proportional scaling factor (see Mollney et 
al. [9] for further details). 

Our Contribution 

This article introduces a new set of variables for simu- 
lating 13 C isotope labeling experiments. The main idea 
underlying this reformulation is that, instead of treating 
fluxes and isotopomer variables separately, we identify a 
set of "isotopically labeled fluxes" as our state variables 
of interest. We refer to these variables as fluxomers. 
Fluxomers combine flux variables with isotopomer vari- 
ables and consequently reduce the complexity and non- 
linearity of the original isotopomer balance equations. In 
this article, we show that by reformulating the flux esti- 
mation problem in terms of fluxomer variables, it is pos- 
sible to construct an algorithm that has the following 
key benefits: 

♦ Provides efficient computation of all isotopomers 
in a metabolic pathway 

♦ Is robust to measurement noise (i.e., suppresses 
the effects of measurement errors) and initial 
conditions 

♦ Eliminates the need for measurement scaling factor 
estimation 

♦ Poses the problem using simple mathematical 
expressions, allowing the use of generic optimization 
algorithms 

The rest of the article is constructed as follows. The 
Results and Discussion section illustrates the advantage 
of our approach via simulation results comparing fluxo- 
mer variables to the commonly used cumomer approach 
and the more recently introduced EMU approach. The 
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Methods section presents the detailed formulation of the 
fluxomers optimization problem and the fluxomers 
iterative algorithm (FIA) that provides a reliable and 
efficient method for solving it. All source code and 
executables for our algorithms are freely available at the 
authors website [16]. 

Results and Discussion 

We compared our FIA algorithm to the widely used 
MFA software 13CFLUX [17], which relies on the 
cumomer approach, and to the more recent Open- 
FLUX [15] software, which is based on the EMU [14] 
approach. In order to compare the methods, we con- 
ducted flux estimations for various well-studied meta- 
bolic pathways. Our first example is based upon the 
tutorial which Wiechert et al. provide with their 
13CFLUX software: the Embden-Meyerhof zxvd Pentose 
Phosphate metabolic pathways of Escherichia coli [17]. 
This example compares the running time and robust- 
ness of both algorithms in response to input noise. Our 
second example compares the results and performance 
of FIA to both an adhoc method and the OpenFLUX 
algorithm for the analysis of lysine production by C. glu- 
tamicum, as described by Becker et al. [18] and Quek et 
al. [15]. 

FIA vs. 13CFLUX Comparison: Embden-Meyerhof and 
Pentose Phosphate Pathways 

In this section we examine a network representing the 
Embden-Meyerhof 'and Pentose Phosphate pathways of E. 
coli, which is based upon the tutorial supplied by Wie- 
chert et al. as part of their 13CFLUX software package. 
Since our FIA implementation natively supports 
13CFLUX input files (i.e. "FTBL" files), the same input 
files can be used for both algorithms. (Note, however, 
that FIA does not require definition of free fluxes nor 
initial values, and thus these are simply ignored when 
imported). Figure 1 shows the simple network used 
along with the nomenclature used in previous publica- 
tions. In addition to the network structure, the models 
are provided with flux and isotopic measurements as 
shown in Table 1. 

First, we examined the output of the two algorithms 
for the traditional "noiseless" input file. In order to run 
the analysis, 13CFLUX requires the user to define a set 
of "free fluxes" along with their associated initial values 
[7]. Note that a bad choice of free fluxes or their asso- 
ciated values can result in poor algorithmic performance 
(both in computation time and accuracy). In fact, under 
various initial guesses the algorithm did not converge at 
all. As for FIA, none of the above is required. Since the 
network along with the given measurements are well 
defined, in the noiseless case the two algorithms 
returned similar values for unidirectional fluxes, as can 




PYR 

Figure 1 E.Coli EMP and PPP Metabolic Pathways. The Embden- 
Meyerhof and Pentose Phosphate metabolic pathways of Escherichia 
coli. 

\ J 



be seen in Table 2. Some slight disagreements were 
observed for the bi-directional fluxes, which are more 
poorly identified. 

We next compared the algorithms' sensitivities to 
noise. In a series of 10 experiments, white Gaussian 

Table 1 EMP & PPP simulation data 



Label input data 



Flux name 


Cumomer Index 


Value 


STD 


GLC 


#000000 


0.445 






#100000 


0.500 






#000001 


0.011 






#000010 


0.011 






#000100 


0.011 






#001000 


0.011 






#010000 


0.011 




Rul5P 


#1xxxx 


0.1979 


0.002 




#x1xxx 


0.0153 


0.002 




#xx1xx 


0.0284 


0.002 




#xxx1x 


0.0122 


0.002 




#xxxx1 


0.0976 


0.002 


Ery4P 


#1xxx 


0.0568 


0.002 




#x1xx 


0.0229 


0.002 




#xx1x 


0.0118 


0.002 




#xxx1 


0.0704 


0.002 


GA3P 


#1xx 


0.0330 


0.002 




#x1x 


0.0126 


0.002 




#xx1 


0.1207 


0.002 


PEP 


#1xx 


0.0330 


0.002 




#x1x 


0.0126 


0.002 




#xx1 


0.1207 


0.002 



Values are taken from the example input file included in the 13CFLUX demo. 
Substrate enrichment values are considered as constants. 
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Table 2 Comparison of FIA with 13CFLUX for the simple 
E.coli metabolic network 



Flux name FIA 13CFLUX 





Est. flux 


MSE 


Est. flux 


MSE 


emp1 


0.5100 


0.0020 


0.5099 


0.0023 


eimp2 


0.8500 


0.0008 


0.8500 


0.0007 


emp3 


0.8500 


0.0008 


0.8500 


0.0007 


emp4 


1 .8700 


0.001 1 


1.8700 


0.0006 


emp5 


1 .8700 


0.001 1 


1.8700 


0.0006 


emp6 


1 .8700 


0.001 1 


1.8700 


0.0006 


ddd1 


0.5100 


0.0019 


0.5101 


0.0023 


ddd2 


4.4234 


0.5483 


4.3281 


0.9652 


ppp2r 


4.0834 


0.5485 


3.9880 


0.9657 


ppp3 


4.4689 


1 .0365 


2.7370 


1.1057 


ppp3r 


4.2989 


1 .0368 


2.5670 


1.1057 


ppp4r 


4.0768 


0.3643 


4.1740 


1.1608 


ppp4 


4.2468 


0.3640 


4.3440 


1.1604 


ppp5r 


0.2538 


0.1535 


0.2680 


0.0654 


ppp5 


0.4238 


0.1531 


0.4381 


0.0655 


ppp6r 


0.2550 


0.0175 


0.2560 


0.0194 


ppp6 


0.4250 


0.0171 


0.4260 


0.0188 


upt 


1 .0200 


0.0004 


1.0200 


0.0001 


coOut 


0.5100 


0.0019 


0.5101 


0.0023 



Comparison of estimated fluxes and mean-square estimation error using 
"noiseless" data. 



noise was added to all of the measured isotopomer 
values, and the outputs and computation times for both 
algorithms were recorded. As can be seen in Figure 2, 
unidirectional fluxes remain quite constant and hardly 
suffer from the added experimental error (for both algo- 
rithms). However, the bi-directional fluxes are affected 
by the added noise. 13CFLUX suffers from a higher var- 
iance spread of the estimated values than FIA (thus is 
more sensitive to the added measurement noise). Note 



Table 3 Algorithm running time comparison for FIA vs. 
13CFLUX 

FIA 6.63 7.56 5.17 6.85 8.83 5.92 9.53 6.47 6.97 6.77 
13CFLUX 59.14 56.93 76 121 65.7 451 81.7 173 177 69.65 
Running time is shown in seconds. 

that the difference arises not only due to the mathemati- 
cal model used, but also due to the stability properties 
of the optimization method chosen. 

We next examined the computational performance of 
the two methods. Table 3 shows the algorithm running 
time for convergence (in seconds). The average running 
time for 13CFLUX was 133 seconds, while for FIA this 
time was 7 seconds. The running time ratio (13CFLUX/ 
FIA) for individual experiments varied between x9 to x75. 

FIA vs. OpenFLUX Comparison: Lysine Production by C. 
glutamicum 

In this section we examine the analysis of the central 
metabolism of two lysine-overproducing strains of Cory- 
nebacterium glutamicum: ATCC 13032 (lysC^ 1 ) and its 
P EFTU fbp mutant. Both express feedback-resistant iso- 
forms of the aspartokinase enzyme lysC, while the latter 
is additionally engineered to overexpress the glycolytic 
enzyme fructose-l,6-bisphosphatase. The example is 
based upon the measurements provided by Becker et al. 
[18], who implemented an ad-hoc program to estimate 
the values of various metabolic fluxes. In their more 
recent article introducing the OpenFLUX software pack- 
age [15], Quek et al. chose to compare their results to 
those of Becker et al. Therefore, we will expand upon 
their comparison using our FIA implementation. The 
input file for FIA was constructed using the measure- 
ments and pathway structure given in [18] and [15]. As 
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ppp2 
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ppp2r 


11 


ppp3 


11 




ppp3r 


11 


ppp4 


11 


ppp4r 


11 


ppp5 


11 


ppp5r 


11 


ppp6 


11 


ppp6r 
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Figure 2 Measured fluxes values. Bidirectional 


fluxes calculated using FIA and 13CFLUX for noisy measurement set. 
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described in [15], the published mass isotopomer frac- 
tions were modified for mass interference from non-car- 
bon backbone isotopes using the molecular formula of 
the amino acid fragments. FIA supports automatic gen- 
eration of the naturally occurring isotopes correction 
matrix when the measured molecular formulas are sup- 
plied. This adjusts the measured fluxomers vector 
appearing in the objective function during the process 
of optimization. If necessary, it is possible not to use 
this feature but instead to directly supply the algorithm 
with the corrected measurement values. 

When comparing the running times of FIA with 
OpenFLUX, the different algorithmic approaches of the 
two must be kept in mind. While OpenFLUX requires 
the user to supply it with sets of free fluxes, FIA 
requires no free fluxes nor initial values. Open-FLUX 
rapidly evaluates dozens of different optimization cycles 
with random initial values and seeks the best fitting 
result among them, while FIA uses only one single 
(longer) run. As such, the convergence probability of 
OpenFLUX depends on the number of attempts and 
random values generated during its operation, while the 
FIA results do not depend on any random value. 
Furthermore, in its analysis, EMU based algorithms eval- 
uate only the fluxes necessary for measurement compar- 
ison, and thus their running time depends both on the 
metabolic network structure and the amount and loca- 
tion of the given measurements. FIA, on the other hand, 
can supply the entire set of metabolic fluxes at any 
given time, with no additional computation requirement 
(which depends mainly on the network structure). 
Measured fluxes as constants 

First, we ran the exact same simulation as Quek et al. 
performed in their article. They supply very accurate 
(mean error in the order of 0.15 mol%) values for the 
label measurements, and used the given measured fluxes 
as if they were noiseless measurements (thus as con- 
stants). We start by comparing the simulation time for 
this simple case. According to [15] and as validated by 
us using our computer, OpenFLUX required 50 itera- 
tions of about 16 seconds each in order to find a decent 
minimal point, hence about 800 seconds in total. While 
so, the FIA analysis took 60 seconds for initial analysis 
and matrices creation, and 300 further seconds for con- 
vergence, thus 360 seconds as a whole. Regarding the 
simulation results, as one can see in Table 4 and Table 
5 the fluxes are very close to those calculated before, 
and the estimated fluxes FIA returned had the lowest 
residual value compared to the other methods. 
Measured fluxes as measurements 

We can also run the same optimization, but weight the 
given flux measurements by their variances. When run- 
ning this optimization using OpenFLUX (again using 50 
iterations), the amount of time was greatly increased, 



and ended in around 48 minutes. For FIA, on the other 
hand, the running time was the same as before, thus 
about 6 minutes. Comparing the results of the algo- 
rithms, OpenFLUX suffered from severe convergence 
problems. Most of its iterations ended without conver- 
ging at all, while those that did converge yielded useless 
results, far from the measurements. FIA, on the other 
hand, succeeded in converging for all scenarios. For the 
wildtype lysine producing pathway, the results were very 
close to the ones before (since the fluxes and measure- 
ments were quite accurate). For the mutant example, 
which was less accurate, a reduction of the residual 
value was achieved by small changes to the measured 
fluxes, fluxes and residual values can be examined in 
Table 4 and Table 5. 
Using non-normalized MS measurements 
We now show that FIA can easily use incomplete or 
non-normalzied measurements by examining its perfor- 
mance in the example above. The supplied MS measure- 
ments were normalized to the n +1 backbone carbon 
atoms of the measured metabolites. Instead of using the 
supplied normalized data, we multiply each set of meta- 
bolite measurements by a random constant number. By 
doing so, we simulate the case in which only the first 3 
(2 for GLY) MS peaks were measured, and had not 
been normalized. The original and supplied non-nor- 
malized measurement values can be found in Table 4. 
Note that the values were corrected by the molecular 
formulas of the measured fragments (again, can be auto- 
matically performed by FIA). In the absence of normal- 
ized data, FIA gave estimated fluxes very close to the 
previous cases, with very low residual values, as can be 
seen in Table 5. The running time of the algorithm was 
not affected by the change. 

Conclusions 

The main contribution of this article is the introduction 
of fluxomers-z. new set of state variables used to simu- 
late 13 C metabolic tracer experiments. The fluxomers 
approach allows the central optimization problem of 
MFA to be reformulated as a sequence of quadratic pro- 
grams, which form the basis of the fluxomers iterative 
algorithm (FIA). Both fluxomers and FIA result in sev- 
eral important benefits compared to flux-isotopomer 
variables. Among these advantages are (i) a reduction in 
algorithm running time required for simulation of isoto- 
pomer distributions and metabolic flux estimation, (ii) 
reduced sensitivity to measurement noise and initial flux 
values and (iii) availability of complete isotopomer infor- 
mation for a given network (as opposed to the EMU 
approach, which only supplies partial information) with- 
out the need for user specification of free fluxes or 
initial flux values. Additionally, the error model used by 
the FIA algorithm has the advantage that it depends 
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Table 4 Relative mass isotopomer fractions comparison for wild-type and mutant C. glutamicum 

Wildtype Mutant 
Fragment Non-normalized Exp. Ad-hoc OpenFLUX FIA Exp. Ad-hoc OpenFLUX FIA 















const. 


meas. 


ratios 








const. 


meas. 


ALA 260 


MO 


206.3562 


0.5085 


0.509 


0.509 


0.5099 


0.5099 


0.5097 


0.5230 


0.525 


0.525 


0.5247 


0.5247 




Ml 


102.8634 


0.3529 


0.354 


0.354 


0.3534 


0.3534 


0.3537 


0.3410 


0.342 


0.342 


0.3425 


0.3425 




M2 


4.8452 


0.1058 


0.106 


0.106 


0.1063 


0.1063 


0.1062 


0.1030 


0.104 


0.104 


0.1037 


0.1037 


VAL 288 


MO 


41.4005 


0.3455 


0.348 


0.348 


0.3459 


0.3458 


0.3457 


0.3640 


0.366 


0.366 


0.3661 


0.3663 




M1 


39.6134 


0.3983 


0.398 


0.398 


0.3986 


0.3986 


0.3987 


0.3920 


0.392 


0.392 


0.3921 


0.3922 




M2 


10.7340 


0.1845 


0.184 


0.184 


0.1846 


0.1846 


0.1847 


0.1750 


0.175 


0.175 


0.1750 


0.1749 


THR 404 


MO 


194.9082 


0.3330 


0.334 


0.334 


0.3343 


0.3343 


0.3340 


0.3440 


0.344 


0.344 


0.3439 


0.3439 




M1 


159.2226 


0.3764 


0.376 


0.376 


0.3757 


0.3757 


0.3759 


0.3730 


0.371 


0.371 


0.3715 


0.3721 




M2 


35.2094 


0.1957 


0.196 


0.196 


0.1956 


0.1956 


0.1957 


0.1910 


0.192 


0.192 


0.1920 


0.1918 


ASP 418 


MO 


159.9111 


0.3343 


0.333 


0.333 


0.3337 


0.3337 


0.3334 


0.3450 


0.343 


0.343 


0.3432 


0.3433 




M1 


128.3755 


0.3732 


0.375 


0.375 


0.3750 


0.3750 


0.3752 


0.3700 


0.370 


0.371 


0.3708 


0.3714 




M2 


28.7782 


0.1955 


0.196 


0.196 


0.1960 


0.1959 


0.1960 


0.1920 


0.193 


0.192 


0.1924 


0.1922 


GLU 432 


MO 


3.8009 


0.2469 


0.25 


0.249 


0.2474 


0.2473 


0.2469 


0.2570 


0.264 


0.264 


0.2634 


0.2624 




M1 


4.4232 


0.3648 


0.366 


0.366 


0.3661 


0.3661 


0.3660 


0.3650 


0.365 


0.365 


0.3656 


0.3658 




M2 


1 .7429 


0.2412 


0.239 


0.240 


0.2406 


0.2406 


0.2409 


0.2360 


0.232 


0.232 


0.2322 


0.2327 


SER 390 


MO 


224.9043 


0.4497 


0.449 


0.448 


0.4487 


0.4488 


0.4490 


0.4620 


0.463 


0.463 


0.4635 


0.4628 




M1 


108.4056 


0.3576 


0.358 


0.358 


0.3578 


0.3578 


0.3580 


0.3490 


0.349 


0.349 


0.3491 


0.3492 




M2 


3.5199 


0.1428 


0.143 


0.144 


0.1437 


0.1437 


0.1434 


0.1400 


0.140 


0.140 


0.1399 


0.1403 


PHE 336 


MO 


250.7079 


0.2712 


0.274 


0.274 


0.2764 


0.2764 


0.2769 


0.2870 


0.289 


0.289 


0.2881 


0.2874 




M1 


303.6304 


0.3816 


0.381 


0.381 


0.3817 


0.3817 


0.3822 


0.3800 


0.381 


0.381 


0.3809 


0.3806 




M2 


129.5861 


0.2282 


0.228 


0.228 


0.2263 


0.2264 


0.2261 


0.2200 


0.220 


0.220 


0.2206 


0.2210 


GLY 246 


MO 


738.7580 


0.7407 


0.742 


0.742 


0.7417 


0.7417 


0.7421 


0.7410 


0.743 


0.743 


0.7426 


0.7426 




M1 


39.7395 


0.1845 


0.185 


0.185 


0.1852 


0.1852 


0.1849 


0.1830 


0.184 


0.184 


0.1844 


0.1844 


TYR 466 


MO 


36.7321 


0.2344 


0.236 


0.236 


0.2380 


0.2380 


0.2384 


0.2460 


0.249 


0.249 


0.2481 


0.2475 




M1 


43.7966 


0.3530 


0.356 


0.356 


0.3567 


0.3567 


0.3572 


0.3510 


0.358 


0.357 


0.3572 


0.3569 




M2 


18.6839 


0.2423 


0.245 


0.245 


0.2433 


0.2433 


0.2431 


0.2340 


0.238 


0.238 


0.2387 


0.2390 


TRE 361 


MO 


34.1048 


0.0613 


0.062 


0.062 


0.0612 


0.0612 


0.0608 


0.0880 


0.088 


0.088 


0.0884 


0.0884 




M1 


327.3441 


0.6040 


0.607 


0.606 


0.6051 


0.6051 


0.6057 


0.5730 


0.577 


0.574 


0.5743 


0.5742 




M2 


27.0318 


0.2070 


0.207 


0.207 


0.2084 


0.2084 


0.2084 


0.2130 


0.213 


0.213 


0.2128 


0.2126 


Sum of weighted residuals 




761 


684 


654 


650 


718 


1735 


1461 


1451 


1308 





Experimental and calculated isotopomer MS fractions. The experimental data and ad-hoc simulation results are taken from Becker et al. [18]. The OpenFLUX 
results are taken from [15]. The simulated "non-normalized" data is generated by multiplying the given values after natural isotope correction by random factors. 
Several FIA estimations are provided: using the given fluxes as constants (under "const"), as measurements (under "meas."), and when using the simulated non- 
normalized data (under "ratios"). As can be seen, FIA agrees with previous results (even when the data is used without normalization). For the mutant case, 
better fits are achieved when allowing the supplied fluxes to change as well. 



solely upon isotopomer ratios rather than complete iso- 
topomer fractions, and therefore it eliminates the need 
to estimate a normalization factor for each measured 
isotopomer distribution. Our current results show signif- 
icant improvements even with regards to simplistic tra- 
cer experiments (the running times have been improved 
by an order of x3 to x20 compared to the 13CFLUX 
algorithm, and about x2 to x8 compared to the Open- 
FLUX implementation). It is important to note that the 
total time required to obtain an MFA solution is con- 
trolled both by (i) the time of each iteration and (ii) the 
number of optimization iterations that are required to 
achieve a reliable solution. While a single OpenFLUX 



iteration is certainly faster than a single iteration of FIA, 
the FIA algorithm was expressly constructed to provide 
high reliability in achieving the optimal solution. There- 
fore, FIA was able to consistently find a better optimal 
solution in less total time in comparison to the other 
algorithms examined. Furthermore, extending the fluxo- 
mers formulation to other global optimization techni- 
ques is straightforward. We expect that reformulating 
more sophisticated MFA problems-for example, invol- 
ving optimal experimental design or large-scale meta- 
bolic networks-in terms of fluxomer variables will lead 
to dramatic enhancements of algorithmic efficiency and 
robustness. 
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Table 5 Metabolic fluxes comparison for wild-type and mutant C. glutamicum 

Wildtype Mutant 



Becker OpenFLUX FIA Becker OpenFLUX FIA 









const. 


meas. 


ratios 






const. 


meas. 


Glucose 6-phosphate isomerase 


49.8 


51 .2 


51 .9 


52.0 


51 .5 


41 .6 


40.4 


42.1 


42.5 


Glucose 6-phosphate dehydrogenase 


46.8 


45.0 


44.7 


44.7 


45.1 


56.2 


57.5 


55.7 


55.1 


TrarKalHnl^KP 

1 1 CM 1 JO 1^ 


14 


13.4 


1 3.3 


13.3 


13.4 


17.5 


1 7.7 


1 7.3 


1 7.0 


Transketolase 1 


14 


13.4 


13.3 


13.3 


13.4 


17.5 


17.7 


17.3 


17.0 


Transketolase 2 


11.9 


11.3 


11.2 


11.2 


11.3 


15.8 


16.4 


15.6 


15.4 


Glyceraldehyde 3-phosphate dehydrogenase 


157.5 


158.0 


158.2 


158.6 


158.0 


160.8 


161.0 


161.0 


160.5 


Pyruvate kinase 


147.3 


148.0 


147.8 


148.2 


147.6 


152.6 


152.0 


152.5 


152.0 


Pyruvate dehydrogenase 


77.5 


75.8 


74.8 


74.9 


74.9 


87.5 


85.2 


85.1 


79.7 


Pyruvate carboxylase - carboxykinase 


34.4 


35.8 


35.9 


36.1 


35.8 


31.5 


32.4 


32.5 


34.9 


Citrate synthase 


52.5 


50.8 


49.6 


49.7 


49.9 


67.7 


65.4 


65.3 


58.9 


Isocitrate dehydrogenase 


52.5 


50.8 


49.6 


49.7 


49.9 


67.7 


65.4 


65.3 


58.9 


Oxoglutarate dehydrogenase 


41.2 


39.4 


38.2 


38.3 


38.5 


59.9 


57.6 


57.5 


50.7 


Aspartokinase 


11.2 


11.2 


11.2 


11.4 


11.2 


14.2 


14.2 


14.2 


15.9 



Estimated metabolic fluxes values for the different approaches - the ad-hoc simulation results from Becker et al. [18], the OpenFLUX results [15], and the FIA 
results for its various simulated scenarios (measured fluxes used as constants, as measurements, and when using ratios of non-normalized data.) 



Methods 

In the following, we show how to construct and solve 
MFA problems using fluxomer variables. First we 
define and explain the basic properties of fluxomers. 
Then we show how to express MFA balance equations 
and measurements in terms of fluxomers. Finally, we 
formulate the MFA optimization problem and present 
the FIA algorithm for solving it. Throughout this sec- 
tion we use boldface uppercase letters A to denote 
matrices, lowercase boldface letters x to denote vec- 
tors, and lowercase letters u for scalars. We use the <° 
>product z = xQy to represent the element-wise pro- 
duct vector, i.e. z t - Xij^ The model formulation will 
be illustrated using the simple metabolic network 
shown in Figure 3. 

Fluxomers overview 

Traditional MFA approaches construct distinct variables 
for each flux and for each possible labeling state (isoto- 
pomer) associated with all metabolites in the network. 
Fluxomers, on the other hand, are a composite of these 
two and therefore allow the network state to be 
described using only one variable type. 

Definition 1 (Fluxomer) A fluxomer is the rate that a 
metabolic reaction transfers labeling from one or more 
specific substrate isotopomers into product isotopomers. 

Taking each fluxomer to be a transformation from 
one set of labeled atoms into others, we can write its 
labeling state as an array of binary elements represent- 
ing the state of each atom it consumes (0 representing 
an unlabeled atom and 1 representing a labeled atom). 
Thus, /(lOOl) is a fluxomer of reaction i consuming 4 



atoms, with its first and last atoms labeled and two mid- 
dle atoms unlabeled. When using x as an index for one 
(or more) of the atoms, we denote a sum of fluxomers 
where the indicated atom can be either labeled or unla- 
beled (e.g., 7X1*01) is the sum of/;(1001) and^(HOl)). 
See Figure 3b for a detailed example. 

Traditional metabolic fluxes and isotopomer variables 
can be easily expressed using fluxomers. We start with 
metabolic fluxes, which are just a sum of their 




(b) 



.upt 




Figure 3 Simple metabolic network, (a) Standard network 
representation. Carbon atoms are drawn explicitly with arrows to 
indicate atom transitions. Unidirectional arrows represent 
unidirectional fluxes while bidirectional fluxes (such as flux 5) are 
represented by bidirectional arrows, (b) Fluxomers representation. 
Each arrow is a group of fluxomers. X's appear on the appropriate 
atom positions to indicate summation of divergent fluxomers. 
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associated fluxomers. For the simple network in Figure 



3b we have: 




fupt =fl = 


y^/i (ijkl) = fi (xxxx) 


h = 


^fi{ij) =fi{xx) 


fs = 




fout = fl = 


y^/4(ijfe/) = f^{xxxx) 


fs = 


^/5(zjfe/) = fs(xXXX) 


fsr = 


^2f5r{iftl) = f 5r {xXXX) 



(1) 



We can also express isotopomer abundances in terms 
of fluxomer variables for the same example. Because of 
the assumption that enzymes do not differentiate 
between the various isotopomers of a given metabolite, 
the isotopomers within each metabolite pool are distrib- 
uted uniformly across the outgoing fluxes emanating 
from that pool Therefore, the fractional abundance of a 
given isotopomer within a metabolite pool will deter- 
mine the fractional contribution of its corresponding 
fluxomers to the fluxes leaving that pool: 



Aijfe/ =f\{i]kl)lh 



Q 
Eijki 



ij -f4{xxij)/f 4 =f 5 r{ij)/f5r 



■■ f 4 {ijxx)/f 4 = 

■■ fout{i)kl)/fout 



mm 
mm 



(2) 



Fluxomer balance equations 

We now examine the fluxomer balance equations that 
describe how fluxomers are propagated through the 
metabolic network. These balance equations represent 
the main mathematical device for calculating steady- 
state fluxomer values for a given network. For ease of 
notation, let us define the vector of metabolic fluxes in 
our system by u £ lZ n and the vector of fluxomers as 
x e 1Z m . As shown above, the metabolic fluxes are calcu- 
lated from a linear transformation of the fluxomers. 
Denoting this linear transformation matrix as U, we can 
write u = Ux. We now assume that we are given a cer- 
tain u vector and wish to calculate the fluxomers in our 
system. We start by considering balances on "simple 
fluxomers", i.e. those that originate exclusively from a 
single metabolite pool. (An example of a simple fluxo- 
mer is ^(Ol) in Figure 3, which derives solely from pool 
B.) Under conditions of metabolic and isotopic steady 
state, the rate of 01-labeled molecules entering pool B 
must balance the rate that 01-labeled molecules leave 
that pool. Therefore, we can construct a balance on 
fluxomers around pool B as 



However, according to eq. 2 the left-hand side of this 
equation can be re-expressed as 

Boiifs + fi) = ^r^ifs + h\ Substituting this latter 

h 

result into the flux balance equation and solving for the 
fluxomer / 5 (01) yields 



/s(01) 



fs 



fs+fi 
2(u)(h T x), 



(fi(01«)+/3(01)+/ 5r (01)) 



(4) 



where g(u) is a function of u alone, and h is a con- 
stant vector. Thus, for this simple case we can solve for 
the outgoing fluxomer /s(01) directly in terms of the 
fluxomers entering pool B and the total fluxes f 2 and f 5 
leaving pool B. 

We now turn to the more complex situation in which 
the output fluxomer originates from more than one 
metabolic pool. For example, consider fluxomer ^(OOOl) 
coming from pools C and D. Here, the fraction of 0001- 
labeling carried by flux f 4 is proportional to the abun- 
dance of 01 -labeling in C and 00-labeling in D: 



/ 4 (0001) = /4C01D00 

/ A(^oi) + / 5 (oi) \ (Moon 

V U /V/3+/4/ 



(5) 



As before, the outgoing fluxomer / 4 (0001) can be 
expressed solely in terms of g-& pure function of u 
(always a rational function of outgoing fluxes)-and a 
product of linear projections of x. 

Without loss of generality, we restrict ourselves to 
fluxes coming from at most 2 metabolic pools (referred 
to subsequently as the "left" and "right" pools). When 
the system reaches steady state, we have 



x = g(u) o (Hix) o (H 2 x), 



(6) 



/ 5 (01) +/ 2 (01) =/i(01xx) +/ 3 (01) +/ 5r (01). 



(3) 



where g is a function TU 1 —> 7Z m , and (H lf H 2 ) are two 
m x m matrices. This equation allows for the output 
fluxomers emanating from a specific metabolite pool to 
be expressed in terms of the total flux vector u and the 
fluxomers entering the pool. This enables each outgoing 
fluxomer to be solved "locally" for the incoming fluxo- 
mers. Note that this local calculation does not involve 
any matrix inversions or other expensive computational 
procedures. If there are no recycle loops in the network 
so that all possible paths through the network are non- 
selfintersecting, this equation can be used to solve 
sequentially for all "downstream" fluxomers in terms of 
previously calculated "upstream" fluxomers. In the pre- 
sence of recycle loops an iterative approach can be con- 
structed to solve for the fluxomers while still avoiding 
repeated matrix inversions. 
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Constructing the system matrices 

The matrices H lf H 2 e JZ mxm are defined by 

(Hi)y = 1, if Xj enters the left (for H 2 , right) 

source metabolic pool in a reaction for which X{ is a product (7) 
(Hi)^- = 0, otherwise. 

The function g : 7Z m -> lZ n is defined as 



ft(u) 



T 

Sii u 



with g!i, g 2i , g 3i e 1Z n given by 
(Sii)j = 1' ^ tne fluxomer Xj is part of the flux/j, 

= 0' otherwise 
(g 2i ) ; = 1, if flux /J exits the left source pool 
in a reaction for which X{ is a product, 
(g 2i ) ; = 0, otherwise 

(g 3i ) ; = 1, if flux/j exits the right source pool 

in a reaction for which x\ is a product, 
(foi)j = 0/ otherwise. 

In matrix form, 

Giu 



(8) 



(9) 



8(u) 



(G 2 u)o(G 3 u)' 



(10) 



Isotopomer measurement formulation 

In the following, we develop a systematic method for 
expressing measured isotopomer variables using fluxo- 
mer notation. The final result of the analysis shows that 
isotopomer measurements can be written simply as the 
norm of a linear transformation of fluxomers, thus Err 
~ ||Ax|| 2 . First, we briefly summarize the available iso- 
topomer measurements provided by Nuclear Magnetic 
Resonance (NMR) and Mass Spectrometry (MS) meth- 
ods. We then discuss the mathematical modeling of 
these measurements using fluxomer variables. 
Available isotopomer measurements 

MFA experiments are typically carried out by (i) intro- 
ducing a labeled substrate into a cell culture at meta- 
bolic steady state, (ii) allowing the system to reach an 
isotopic steady state, and (hi) measuring isotopomer 
abundances of metabolic intermediates and byproducts 
using either MS or NMR analysis. These two measure- 
ment techniques provide qualitatively different informa- 
tion about isotopic labeling. 

♦ X H NMR: Measures the fractional 13 C enrichment 
of each proton-bound carbon atom, irrespective of 
the labeling of its neighboring carbon atoms. Both 



12 C and 13 C atoms are distinguishable in the same 
spectrum, and therefore the peak areas corresponding 
to different carbon isotopes can be normalized directly. 

♦ 13 C NMR: Quantifies isotopomers based on the 
presence of multiplet peaks (e.g., doublets, triplets, 
doublet doublets, etc.) in the spectrum caused by 
two or more neighboring 13 C atoms. Because 12 C 
atoms are undetectable by 13 C NMR, it is impossible 
to quantify the overall fraction of each isotopomer 
unless X H NMR spectra are simultaneously obtained. 
Instead, only the relative ratio of different isotopo- 
mers can be assessed by 13 C NMR. 

♦ MS: This technique is usually preceded by some 
form of chromatographic separation (GC or LC) to 
resolve mixtures into their individual components. 
These components are then ionized and fragmented 
in the MS ion source. The ionized particles are sepa- 
rated according to their masses by an electromag- 
netic filter, and a detector measures the relative 
abundance of each mass isotopomer. These abun- 
dances can be normalized to a fractional scale if all 
MS peaks corresponding to a particular ion are 
simultaneously measured. 

Previous studies based on flux-isotopomer variables 
have modeled the measurement error as Gaussian noise 
added to the fractional isotopomer enrichments. There- 
fore, if y is the vector of measured isotopomer fractions, 
this model states that y = y + e, where e is the Gaussian 
error term. However, a more accurate error model 
would add the measurement noise directly to the physi- 
cally measured values. The motivation for the tradition- 
ally chosen error model is its relative simplicity when 
expressed using flux-isotopomer variables. Furthermore, 
since some isotopomers of a specific metabolite may be 
unmeasurable, the isotopomer fractions cannot be 
experimentally determined in many cases. This implies 
the need for an alternative error model that avoids these 
shortcomings. 
Measurement Error Model 

We denote the measured isotopomer abundances by a 
vector rh. For NMR analysis, the elements of rh are pro- 
portional to the areas under the different spectral peaks. 
For MS, they are proportional to the integrated ion 
counts associated with each mass isotopomer. Since rh 
is the measured quantity, the correct error model is an 
addition of Gaussian noise so that rh = m + e> where m 
is the "true" measurement value. The measured isotopo- 
mer fractions y are then expressed as 



Yi 



J2i rhi Ei ( m i + e i) 



(11) 
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Let Sj represent the residual between the modelpre- 
dieted and experimentally measured abundance of a sin- 
gle isotopomer. After multiplying eq. 11 by J2i ( m * + e i) 
and rearranging, the residual expression becomes 

£ i = m i ~ UJ2 m ') = e i - » 12 ei > (12) 

i i 

where Sj is a sum of Gaussian variables. Noting that 
each measurement m ; is simply proportional to a linear 
combination of fluxomers, the residual expression eq. 12 
takes the form 

£ = [diag(y)T-V]x, (13) 

where T and V are transformation matrices needed to 
convert fluxomers to isotopomer measurements and the 
diag operator converts its vector argument into a diago- 
nal matrix. The resulting expression is both a simple 
sum of Gaussian vectors and affine in x. 

The advantage of this objective function is that it only 
depends upon the relative isotopomer intensities in the 
vector y but does not depend upon how these intensi- 
ties have been normalized (as long as the transformation 
matrix T is constructed accordingly). This eliminates the 
need to estimate optimal normalization factors that are 
required by previous algorithms in order to convert 
experimental measurements into isotopomer fractions. 
This is true for both MS and NMR measurements, 
either when conducted alone or used together in the 
same experiment. 

The MFA optimization problem using fluxomers 

Now that we have defined both the isotopomer mea- 
surements and the feasible solution set, we can formu- 
late the least-squares MFA optimization problem in 
terms of fluxomer variables. Our objective is to find the 
flux vector u that minimizes the measurement error. In 
addition to the fluxomer balances, usually upper bounds 
u ub are provided for all fluxes. As has been proven by 
Wiechert et al. [6-9], once the inputs to the system and 
u are set, the solution (x, u) is unique. In other words, 
the steady-state fluxomer balance equation, eq. 6, is 
actually an implicit definition of x(u). With this in 
mind, the MFA optimization problem can be simply 
defined as 

min||Ax(u) -b|| 2 (14) 

ugQ 

with 



where A selects the measured elements of the fluxo- 
mers vector x(u), b contains their associated values, and 



S is the stoichiometric matrix of the reaction network. 
Note that b may contain non-zero elements only when 
associated with measurements of absolute flux values. 
For isotopomer measurements, the associated elements 
of b are zero. 

Eq. 14 can be solved using various non-convex global 
optimizing techniques. These optimizers typically 
require the user to provide subroutines for computing 
the value of the objective function and its first deriva- 
tives at various points along their convergence path. 
Furthermore, evaluation of the function x(u) and its 
derivatives are the main (practically only) time-consum- 
ing procedures when solving the optimization problem 
in eq. 14. The mathematical formulation of eq. 14 is 
similar to the optimization problem resulting when 
using the labels and fluxes variables, with one exception 
- the implicit formula for x(u). As shown above, using 
fluxomers we are able to formulate the propagation 
equation (and thus solving x(u)) as a multiplication of 
homogeneous functions of fluxes, and second order 
functions of fluxomers. Using labels and fluxes, formu- 
lating the same equation results in a sum of functions of 
the same structure, and the homogeneous separation 
property vanishes. The following sections exploit this 
unique property of the fluxomers propagation equation 
in order to achieve great reduction in the system com- 
putational complexity, leading to the FIA algorithm. 

Fluxomers Iterative Algorithm (FIA) 

This section deals with the evaluation of x(u) along with 
its gradient using the fluxomer formulation. First, we 
show that x(u) can be calculated iteratively while avoid- 
ing repeated matrix inversions. Then, we demonstrate 
how the number of iterations can be reduced using a 
Newton-type gradient-based algorithm. Finally, we 
explain how it is possible to greatly increase the sparsity 
of the system using a simple linear transformation of 
variables, which further reduces the number of iterations 
needed for convergence. 
Solving the fluxomer balance equations 
A simple approach for computing x given u is to imi- 
tate nature. Once a metabolic network reaches steady 
state (namely, when u is constant), changing its input 
labeling does not affect its flux values u, but only influ- 
ences the labeling of its intermediate metabolite pools. 
The metabolite labeling patterns become gradually 
mixed and propagated throughout the network until 
isotopic equilibrium is reached. Accordingly, a simple 
approach for solving eq. 6 is by using its iterative repre- 
sentation (which is similar to the process taking place 
in nature): 

x t+ i = g(u) o (Hix t ) o (H 2 x t ), 
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where x t is the fluxomer vector at iteration t and x t 
+1 is the fluxomer vector at iteration t + 1. In order to 
simulate the steady-state labeling, we initialize the sys- 
tem with the vector x 0 in which only the input fluxo- 
mers are labeled and all others are unlabeled. By 
recursively substituting x back into the equation, 
steady state is eventually reached and the final value of 
x is obtained. (This equation represents a non-linear 
time-invariant Markov chain.) For the Emb den- Meyer - 
hof and Pentose Phosphate Pathway example in the 
Results and Discussion section, it takes a few hundred 
iterations to achieve complete stability of the solution 
(maximal fluxomer value change on the order of le- 
17). Algorithm convergence for a given input vector is 
retrieved exactly as in the real biological system, and 
thus a unique solution always exists (for realistic meta- 
bolic networks). 

We now show it is possible to reach pathway conver- 
gence in much fewer iterations. First, we write eq. 6 as 



F(x,u) = g(u) o (Hix) o (H 2 x) - x. 



(15) 



Now, in order to find the values of (x, u) one needs to 
solve F(x, u) = 0 while holding u constant. We choose 
to use one of the classic and powerful algorithms for 
finding roots of an equation, the well known Newton- 
Raphson [19-21] method. Roughly speaking, the change 
of the x vector at each iteration is calculated by 

x t+ i = x t - (Fx(x,u)) _1 F(x,u), 

with F' fx,ul = dF( x ' u ) The main concern now is the 

evaluation of the expression (F x (x,u)) _1 F(x,u} Here, it 
turns out that due to the decomposable nature of F(x, 
u), the derivative F x at a point (x, u) is the simple matrix 



F;(x,u) = (g(u)o(H lX ))H 2 

+ (g(u) o (H 2 x)) Hx - I. 



(16) 



Therefore, finding r = (F^(x,u)) 1 F(x / u) is equivalent 
to solving the linear system of equations 



(F;(x,u))r = F(x,u). 



(17) 



In order to determine the root of the propagation 
equation, FIA starts with an iteration or two using New- 
ton's correction and then continues with the simple 
"natural" approach. Applying this method to the Emb- 
den-Meyerhof and Pentose Phosphate Pathway example 
in the Results and Discussion section, only a few dozen 
iterations are now needed. In the next section we show 
how to reduce both the number of variables and the 
number of iterations required for convergence by 
another order of magnitude, without affecting system 
convergence stability. 



Reducing system complexity 

The following section introduces a mathematical 
approach for reducing the number of nonzero elements 
in our system. Variable reduction techniques such as the 
recently developed Elementary Metabolite Unit (EMU) 
network decomposition [14] were developed for applica- 
tion to systems that are modeled using flux-isotopomer 
variables. Fluxomers and the FIA algorithm, as opposed 
to prior approaches, allow us to effectively reduce the 
number of system variables using a simple linear trans- 
formation on x. Our main goal here is to find a trans- 
formation for the fluxomer vector x, y = Kx that: 

♦ Reduces the number of its nonzero elements. 

♦ Reduces the computational complexity of solving 
eq. 16. 

♦ Eases the evaluation of eq. 15. 

From eq. 16 we see that the greatest expense is due to 
inversion of a sum of two linear transformations (Hx 
and H 2 ) of x. From the iterative propagation equation, 
eq. 15, we see that x is iteratively calculated by comput- 
ing the product of the same two matrices. Had it been 
possible to find a sparse, close-to-diagonal representa- 
tion for both H x and H 2 by simply multiplying them by 
the matrix from the right, both problems would be 
solved. 

In order to acomplish the above, we examine the 
properties of the concatenation of these two matrices 
which we denote by H. Next we find the LU factoriza- 
tion of H, 



H = L H U H 



Lhi 

LH2 



U h , 



(18) 



with L H lower triangular and U H upper triangular 
matrices. The matrix L H1 contains the first m rows of 
L H and L H 2 contains the last m rows of L H . Our new 
set of variables now becomes y = U H x, and the new 
propagation equation is 



u h [g(u) o (L m y) o (L H2 y)] - y = 0. 



(19) 



When expressed in terms of the variable y, our system 
becomes much more sparse. This is illustrated in Figure 4 
which shows H lt H 2 , L H1 , L H2 and U H for the Embden- 
Meyerhof and Pentose Phosphate Pathway example. The 
transformation has two essential benefits: 

1. Reduced computational complexity - note that 
the derivative (F x (x,u)) _1 F(x,u) depends upon the 
matrices H x and H 2 which have already been fac- 
tored, and thus solving Newton's step is 
straightforward. 

2. Fewer iterations needed for convergence. 
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ti 




1000 2000 3000 4000 5000 

nz = 161608 



1000 2000 3000 4000 5000 

nz = 1031 168 



m 

li 



7 



1000 2000 3000 4000 5000 

nz = 8616 



1000 2000 3000 4000 5000 

nz = 7234 



1000 2000 3000 4000 5000 

nz = 10824 



Figure 4 System matrices complexity reduction. H 1; H 2 , L m , L H2 and U H for the simple E. coli example. A substantial reduction in nonzero 
elements between the H and L matrices can clearly be seen. 



As a matter of fact, this transformation reduced the 
number of iterations needed for convergence of the sim- 
ple E. coli example to a total of 5. 

Finding 

As discussed above, our optimization problem seeks 
the minimum of ||Ax(u) - b|| 2 . In order to converge 
rapidly, the gradient of the objective function must be 
computed at each iteration of the algorithm. The key 
step for computing it is the evaluation of |* (the deriva- 
tive of the fluxomers as a function of the metabolic 
fluxes). Since we have an implicit function F(x, u) along 
with a valid solution for F(x, u) = 0, the implicit func- 
tion theorem [22,23] can be used to compute 
Because F(x, u) is continuously differentiate around its 
root, we can write 



ax 



aF^un -1 3F(x,u) 



dx 



du 



(20) 



In the previous section we showed that dF( x ' u ) can ^ e 

dx 

directly expressed in terms of the system matrices and 
the vectors x and u. The same procedure can be carried 

aF(x,u) 



out in order to determine 



3F(x,u) 
du 



diag{{Y\\x) o (H 2 x)) 



(G 2 u) o (G 3 u) 



((du) o (G 2 u))G 3 + ((Gxu) o (G 3 u))G 2 
((G 2 u) o (G 3 u)) 2 



(21) 



Keeping in mind that F x (x, u) is in its reduced form 
due to our variable transformation, solving the equation 
3F(x,u)^ /^9x\ 9F(x,u) 
3u / 3u 



can be accomplished 



3x 

efficiently. 
The initial point 

The generation of the initial point for the FIA algorithm 
is very similar to the standard method used by many 
iterior point algorithms for finding a valid initial point 
over a convex linear set. We added a fluxes-measure- 
ment regularization factor A in order to generate an 
initial point closer to the final estimation (and thus 
speed up the convergence process). The initial point is 
generated by solving the following simple convex opti- 
mization problem: 



min(— s + X 

u,seS 



Au 



•u" 2 



du 
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with 



£ 



Su = 0 

[I,-I]U+[I,I] [0,U b ] T >5 



with A a matrix that selects the measurable elements 
of u, u the meaured elements of u (if they exist), 0 a 
vector of zeros, and u b a vector of the flux upper 
bounds. The regularization factor A starts with some 
large value, and if necessary is reduced until the optimal 
value of s is greater than 0. Note that when X — > 0 the 
problem reduces to finding a feasibile solution of u, and 
thus always has a solution (for well-structured 
networks). 
The algorithm 

Summarizing the above discussion leads to the following 
algorithm for efficient solution of the MFA optimization 
problem using fluxomers: 

i. Matrix preparation (run once per network): 
Hi \ / Lhi 



0. Calculate 
factorization. 



H 2 



-H2 



U H using LU (PQ) 



ii. Call the optimizer in order to solve 

min|| Ax(u)-b|| 2 , Q=(u: ^ U = ° 

ugQ v 7 I 0 < u < u ub 



When requested by the optimizer, return x(u) and 
its first derivative: 

1. Set y 0 = Vinit- 

2. Set i = 1. 

3. Calculate 

Yi = U H [g(u) o (L H iVi „ J o (L H2 yi - i)]. 

4. If Wyt-yt _ i|| 2 >s N 

(a) Solve (F;(x,u))r = F(x,u) 

(b) Set Yi = Yi - r according to Newton's method. 

5. If Wyt-yt - i|| 2 > s e go to 3. 

6. Calculate x = [g(u) o (L HlYi ) o (L m ydl 
9F(x,u)^ fdx\ _ 3F(x,u) 



7. Solve 



ax 



du 



The supplied software uses a variant of the "sequen- 
tial least-squares" algorithm [24,25] for solving the 
non-convex optimization problem in eq. 14. This 
essentially breaks the problem into a sequence of con- 
vex optimization problems for which the solution can 
be readily computed. Note that other algorithms can 
be easily used with the same procedures described 
above. 
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